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We revisit here the stabihty of a deformable interface that separates a fuhy- 
developed turbulent gas flow from a thin layer of laminar liquid. Although this 
problem has been investigated in many previous studies, a model that requires no 
parameters that must be chosen a posteriori, and that uses a base state profile that 
is validated against experiments is as yet unavailable. Furthermore, the significance 
of wave-induced perturbations in turbulent stresses remains unclear. Much emphasis 
in the oceanographic literature has been on the distortion of turbulence by the pres- 
ence of interfacial waves, but the relevance of this has not been tested for turbulent 
flow over a thin liquid layer in industrial channel or pipe flows. Therefore, unlike 
previous work, the turbulent base state velocity profile proposed here requires only a 
specification of a fiowrate or pressure drop, and no a posteriori choice of parameters. 
Moreover, the base state contains sufficient detail such that it allows for instability 
due to a viscosity-contrast mechanism (which turns out to be dominant) as well as 
instability due to a critical-layer-type mechanism, and it is validated against the ex- 
perimental and numerical data available in the literature. Furthermore, the effect of 
perturbations in the turbulent stress distributions is investigated, and demonstrated, 
for the first time, to be small for cases wherein the liquid layer is thin. The detailed 
modelling of the liquid layer elicits two unstable modes, and mode competition can 
occur, although in most cases the instability is due to the viscosity-contrast mech- 
anism. In particular, there is the possibility that surface roughness can reduce the 
growth rate of the interfacial mode, and promote a liquid-layer instability to the 
status of most dangerous mode. Our base-state model facilitates a new definition 
of 'slow' and 'fast' waves. We use our linear stability analysis to determine the fac- 
tors that affect the wave speed and demonstrate that the waves are 'slow' according 
to the definition proposed here. Finally, we compare our results with experimental 
data, and good agreement is obtained. 



INTRODUCTION 



A linear stability analysis of small-amplitude waves on an otherwise flat liquid film would 
provide a powerful tool in understanding and modelling the onset of droplet entrainment 
from a liquid layer by a shearing superposed turbulent gas flow, which has numerous indus- 



trial applications (e.g. Hall-Taylor and Hewitt (1970)). Furthermore, this would serve as 
a benchmark for direct numerical simulations of two-layer flows, as in Boeck et al. (2007); 



Fuster et al. (2009); Valluri et al. (2008, 2010) for laminar flows. Although early work in the 



* Present address: School of Mathematical Sciences, University College Dublin, Belfield, Dublin 4, Ireland. 



2 



area focussed on a 'divide-and- attack' approach (Benjamin, 1958; Miles, 1962), wherein first 
the perturbation in the shear stress exerted on a wavy film would be estimated from mea- 
surements and models in turbulent flow over a wavy solid wall, which would then feed into a 
stability analysis of the liquid layer, the advance of computational methods and facilities has 
enabled one to solve stability problems wherein perturbations in the gas and liquid are fully 
coupled. Various studies have pursued this for liquid films sheared by turbulent flow of a 
gas (Kuru et al. , 1995 Miesen and Boersma, 1995). However, several difficulties have arisen 



in these studies, which prevent such previous work to be of direct use as benchmark tests 
for direct numerical simulations of turbulent stratified channel fiows and industrial appli- 
cations. First, a robust model for a base-state velocity profile that has been tested against 
experiments and numerical simulations is not available, and its detailed modelling turns 
out to be important. Next, the base state should not require specification of any ad- hoc 
parameters or parameters that must be chosen a posteriori: merely the fiow rate or imposed 
pressure drop (along with the physical and geometrical properties) should suffice. Previous 
models lack at least one of these aspects; the base state model proposed here satisfies all of 
these criteria. A second motivation for this study is to ascertain the role of perturbations 
in turbulent stresses, which are caused by the presence of waves. These could feed back to 
the growth rate and speed of the wave even in a linear analysis. A further objective of this 
study is to confirm the physical mechanism that leads to instability, as various mechanisms 
have been proposed for general two-layer fiows, and to relate this to a classification of slow 
and fast waves. In many cases, previous work was conducted in the context of wind-driven 
waves in oceanographic applications, and it remains unclear whether, for instance, a Miles- 
type instability could be approached in stratified channel fiows. In addressing this issue, 
the present work leads to a new, naturally classification of slow and fast waves. We briefly 
expand on these issues here. 

In previous work, a 'lin-log' base-state profile has been used in a boundary-layer setting, 
and the friction velocity f/* was guessed (Miesen and Boersma, 1995). Others have adopted 

Here, we derive 



1995). 



an empirical profile, the validity of which is unclear (Kuru et al. 
a base-state model that contains no free parameters: the friction velocity is determined 
as a function of Reynolds number. We develop our base state starting from a rigorously 



vahdated model (that of Biberg (2007)) 
near-interfacial zone 



and generalize it in order to take account of the 



a significant region of the fiow for determining the viscosity-contrast 
instability. We compare the resulting model with direct numerical simulations and experi- 
ments. Furthermore, our model contains no logarithmic singularities: full modelling of the 
viscous sublayers is provided, in contrast to the model of Biberg (2007). The base state (i.e., 
the fiat interface profiles) is presented in Sec. [Hj 

At least four mechanisms have been reported that relate to the instability of a laminar 
liquid layer sheared by an external turbulent gas fiow. The first kind of instability was iden- 
tified by Miles (1957), and is called the critical-layer instability. The transfer of energy from 
the mean fiow into the wave perturbations is governed by the sign of the second derivative of 
the base-state fiow at the critical layer - the height where the wave speed and the base-state 
velocity match. Another kind of instability was identified by Yih (1967), and is called the 
viscosity- contrast mechanism. Here the instability arises due to the jump in the viscosity 
across the interface. In addition, instability can arise due to direct forcing by turbulent 



pressure oscillations in the gas ( 


Phillips 


1957 


and Miesen 


1996 


Miesen and Boersma, 


1995), 



The so-called internal mode (Boomkamp 



is observed when the bottom layer is lami- 
nar. This mode derives its energy both from the interface, and from conditions in the bulk 
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of the bottom layer. Thus, this mode persists even when the upper layer is void (Miesen 



and Boersma, 1995). It has the characteristics of a Tollmien-Schlichting wave, and depends 



sensitively on the viscosity c ontrast. This sensitivitv is a mplified when the bot tom layer has 



non-Newtonian rheology (6 Naraigh and Spelt ( |2010[ )). Ozgen et al. ( |l998[ ) studied such 



a scenario for power-law fluids, where they find that the internal mode dominates over the 
interfacial (viscosity-contrast) mode, a reversal of the situation observed by Miesen and co- 



workers (Boomkamp and Miesen 1996 Boomkamp et al. 1997; Miesen and Boersma 1995) 



for the Newtonian case. This effect is especially visible at high Reynolds numbers. This 



mode competition is also possible for Newtonian fluids: Yecko et al. (2002) have observed 



it for two-phase mixing layers. Benjamin (1958) described an instability mechanism called 
non-separated sheltering, wherein the viscous stresses near the interface give rise to a faster 
velocity downstream of the wave crest, and thus cause a pressure asymmetry. These effects 



are also observed in the viscosity-contrast mechanism (Boomkamp et al. , 1997), of which the 



Benjamin description can be regarded as a limiting case. Belcher and Hunt (1993) proposed 
that the turbulent stresses could produce the same effect, although in this paper, we find 
that the viscous stresses are dominant in the creation of the instability, at least under the 
thin-film parameter regime studied here. Finally, we note here that approximate mecha- 
nisms such as those identified in Kelvin-Helmholtz-type theories, although perhaps relevant 
in large-amplitude waves, do not arise in the full linearized problem wherein viscous effects 
are fully accounted for (Boomkamp and Miesen 1996). 

In order to assess which of these mechanisms dominates in sheared liquid films by tur- 
bulent gas flow, it is useful to note that the critical- layer waves are typically fast, and the 
viscosity-contrast waves are typically slow (relative to an appropriate scale); thus, we pro- 
pose in this paper to classify waves as fast or slow depending on the kind of instability 
mechanism at work. This is not only a semantic choice: our definition leads to a rule-of- 
thumb for deciding when waves are slow or fast, a rule wherein the wave speed is measured 
relative to the interfacial gas friction-velocity, which in turn is based on the pressure gra- 
dient in the channel. For fast waves, the liquid layer has been treated like a moving wavy 
wall, while for slow waves, detailed understanding of the liquid layer and interfacial viscous 
sublayer layer is required. 



Boomkamp and Miesen (1996) classify interfacial instabilities according to an energy 
budget: an energy law is associated to the dynamical equations, and positive inputs of 
energy are identified with various physical mechanisms of instability. Thus, they place the 
critical-layer and viscosity-contrast instability into a rigorous framework. They comment on 
the equivalence between slow waves and the viscosity-contrast instability, and the equivalence 
between fast waves and the critical-layer instability; we take this analysis further by deriving 
a prediction for the character of the wave based on the pressure gradient, and several other 



base-state parameters. The classification is presented in Sec. II B 



In order to explain a further motivation for the present study, we first summarize the 
averaging approach presumed herein. Consider a large ensemble of realizations of a (three- 
dimensional) pressure-driven turbulent channel flow. The velocity field contains perturba- 
tions due to turbulence, and due to the presence of small-amplitude waves. At any time, 
a Fourier decomposition can be taken of the interface height and, simultaneously, of the 
velocity and pressure fields. These Fourier-decomposed velocity and pressure fields can be 
averaged over the ensemble of realizations (as well as over the spanwise direction). These 
ensemble-averaged velocity and pressure fields are not uni-directional, but are distorted due 
to the presence of the corresponding (normal mode) interfacial wave. Example fields that 
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have been obtained in conceptually the same manner from DNS (albeit for turbulent flow 
over a wavy wall) can be found in the paper of Sullivan et al. ( 2000[ ). In the present study, 
results are presented (and compared) from several RANS models. Now, a significant issue 
here is that such wave-induced perturbation stresses may, in principle, affect the growth rate 
and speed of waves. Questions concerning the importance of these stresses have been much 



debated in the literature (Belcher and Hunt 



et al. 


1995 


Miles 


1957 



, 1993 


Belcher et al. 


1994 


Janssen 


2004 


Kuru 



Previous studies on sheared liquid films have not accounted for 
these effects, hence the significance of these stresses are not known at present. The present 
study does at last provide convincing evidence that these are indeed not important in the 
determination of wave speed and growth rate for sheared thin films. The linear stability 
analysis and assessment of the significance of perturbation turbulent stresses (PTS) are the 
subject of Sec. 



Ill 



In Sec. IV 



we carry out a detailed linear stability analysis to investigate the character 
of the interfacial instability, and to find out which parameters affect the wave speed. We 
also study there mode competition, which can be achieved by devising a roughened inter- 
face, this being the averaged effect of instantaneous pressure fluctuations there due to the 
Phillips mechanism. Finally, we compare the model predictions with experimental data and 
simplified (Kelvin-Helmholtz-based) theories in Sec. |V| 



II. THE FLAT-INTERFACE MODEL AND ITS PROPERTIES 



In this section we derive a base state appropriate for a two-layer system in a channel, 
described schematically in Fig. [Tj This is an equilibrium state of the system, in the sense that 
the mean velocities are independent of time, and the mean interfacial height that demarcates 
the phases is flat. The bottom layer is a thin, laminar, liquid layer, while the top layer is 
gaseous, turbulent and fully-developed. A pressure gradient is applied along the channel. 
The mean profile of the system is a uni-directional flow in the horizontal, x-direction. In the 
gas layer, near the gas-liquid interface and the gas-wall boundary, the flow profile is linear. 



Yaglom 


1971 


Pope, 


2000 


)■ 


profile ( 


Monin and Yaglom 



and the viscous scale exceeds the characteristic length scale of the turbulence (Monin and 

In the bulk of the gas region, the flow possesses a logarithmic 
1971; Pope, 2000). We assume that the gas-liquid interface 



is smooth, although we do take account of surface roughness in Sec. |II C[ and again in 
Sec. UVCl 

The growth rate of the wave amplitude depends sensitively on the choice of mean flow. 
Therefore, it is necessary to derive a mean flow-profile that incorporates the characteristics 
of the flow observed in experiments. In this section, we generalize the model of |Biberg 



(2007) and accurately model the viscous sublayers found in two-phase turbulent flows. This 
is a non-trivial improvement, since accurate descriptions of near- interfacial conditions are 
important in understanding the stability of the interface. The functional form of the derived 
velocity profile enables us to express the wall and interfacial shear stresses as functions of the 
mean pressure gradient. This treatment also enables us to write down a turbulent closure 
scheme using an eddy viscosity, which is constituted as a simple function of the vertical 
coordinate, z. Using the above-mentioned assumptions and approximations, we derive the 
mean flow in each layer. 
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FIG. 1: A schematic diagram of the base flow. The hquid layer is laminar, while the gas layer 
exhibits fully-developed turbulence, described here by a Reynolds-averaged velocity profile. A 
pressure gradient in the x-direction drives the flow. 

A. The mean flow 



The liquid film: The Reynolds-averaged Navier-Stokes (RANS) equation appropriate for 
the liquid film is the following: 

d'^Uh _ ^ ^ Q 
dz'^ dx 

where Ul and denote the liquid mean flow velocity and viscosity, respectively. Although 
the pressure gradient dp/dx is written using partial derivatives, the assumption of hydro- 
static balance means that this pressure gradient is merely a function of x. For, the hydro- 
static balance assumption amounts to the equation —dpj/dz = pjg, in which j = L,G labels 
the phase. Integrating once gives pj = —pjgz + pj (x), hence dpj/dx = dpj/dx. Since this 
derivative is independent of the 2;-coordinate, it is phase-independent too, and we therefore 
drop the phase-label from this pressure gradient in the rest of the work. We integrate Eq. ([T]) 
and apply the following boundary conditions, which correspond, respectively, to continuity 
of tangential stress at the interface and no-slip at the channel bottom wall: 



dUr 



dz 



T-i = PgU!„ Ul i-dL) = 0, 



(2) 



where —di < ^ < is the domain of the liquid fllm. The quantity is the interfacial stress 
and is the interfacial friction velocity on the gas side. We are going to close the model 
by flnding an expression for this friction velocity in terms of the pressure gradient dp/dx. 
The boundary conditions ^ yield the following relation for the mean-flow velocity proflle 
in the fllm, 



Nondimensionalizing on the scale Uq, where 



PL 



PgUq = h 



' z + di 



dp 



-dL<z< 0. 



(3) 



dx 



6 



gives 



Ur 



Rel 



Ken 



6 



h ' 



where the tildes denote dimensionless quantities: U = U /Uq, and z = z/h, and where 



Ren 



-Re* 



Pg 



Pg 



are the Reynolds numbers. The definition of Re^ differs from the definition of the Reynolds 
number in single-phase channel flow, Rep = pGh^\dp/dx\/2iiQ. They are related, however, 
by the formula Re^ = y/2Rep. Note furthermore that 



t7(0) 



(4) 



The gas layer: The RANS equation in the gas is 



Pg- 



dUr 



dz 



+ ttss 



dp 
n + -^z, 
ox 



(5) 



where ttss = —p{u'w') is the turbulent shear stress due to the averaged effect of the tur- 
bulent fluctuating velocities. In channel flows, it is appropriate to model this term using an 
eddy- viscosity model (Monin and Yaglom, 1971[). In mixing- length theory, the eddy viscosity 



depends on the local rate of strain (Bradshaw, 1974), which means that the turbulent shear 



stress depends on the square of the rate of strain. Instead of this standard mixing-length 
theory, we introduce an interpolation function for the eddy viscosity, which mimics the or- 
dinary mixing-length theory near the interface and near the wall, and transitions smoothly 
from having a positive slope near the interface, to having a negative slope near the wall. 
Thus, the turbulent shear stress is linear in the rate of strain, and 



'T'tss — Pt- 



dU, 



G 



dz 



Pt = npchU^^G {z) ipi (z) (1 - z) 



(6) 



where px is the eddy viscosity, f/*„ is the friction velocity at the upper wall z = 1, and where 
G (z), ipi (z), and ip^ (1 — z) are functions to be determined. Here ipi and ipw are interface 
and wall functions respectively, which damp the effects of turbulence to zero rapidly near 
the interface and the wall, while G is an interpolation function designed to reproduce the 
law of the wall near the interface and the upper wall. This interpolation scheme is based on 



the work of Biberg (2007). The precise choice of G and the wall functions is given below. 



choices that are conflrmed by the agreement between our predictions of the base state and 
experiments and numerical simulation. Substituting Eq. ^ into Eq. ^ gives the formula 



Ug{z) 



Ug (0) + T,h 

Ug (0) + nh 



z/h 



1 



hdp, 
Ti dx ' 



s ] ds 



Pg + KpGhU^wG (s) ipi (s) V'w (1 - s) ' 



(7) 
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where R = ti/t^. Non-dimensionalizing and using Eq. Q, this IS 

Ug{~z 



Reo J 



(8) 



The ratio R can be obtained in closed form as follows. Since 



dp 

-T-w + ^ ~ ^) 



these formulas can be equated to give 



dp 

dx 



or, 



1 + 



dp h 
dx Ti 



RCq 

Re^ 



hence, 



\R\ 



1 



-Re, 



We use the following form for the G-function. This function is designed to reproduce the 
logarithmic profile near the interface and near the upper wall, a result that we demonstrate 
below (P. [9|; further evidence of the correctness of this choice is provided when we compare 



(9) 



our predictions for the base-state velocity with experiments (Sec. II D). Thus, 

+ \R\'/' (1 - sf 



G (s) = s{l- s) 



i?2 (1 _ s f + Rs{l-S) + S2 



< s < 1. 



We use a Van-Driest type of formalism (|Pope[ 2000|) for the wall functions ipi and ip^: 



^w(l 



(10) 



where n, Ai, and A„ are input parameters. We are now in a position to determine Re^: it 
is obtained as the zero of the function U (1; Re^,) = 0, or 



/^G / 1 ^ Rei 



\Reo5^ + —^5 + 



-Reo 



Rel 



l-^s]ds 



Reo Jo l + ^G{s)^,{s)i^^{l-s) 



In summary, we have the following velocity profile in the base state: 



IJ-L 



UTz) 



iReo{S'-6') + ^Jz + 6) 



Ret rz 









('-if-)- 







-5 < 5 < 0, 



0< z<l. 



0. 

(11) 

(12) 
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FIG. 2: Characteristics of the base profile for fixed parameter values {^l/ fJ'G-, Pl/ PGidi/h) := 
(m,r, 5) = (100, 1000,0.1), and Rcq = 1000. (a) The mean velocity profile; (6) The mean velocity 
profile in wall units, showing the logarithmic and viscous layers (the viscous layer has a wall-unit 
thickness of approximately 5); (c) The Reynolds stress profile corresponding to the basic velocity. 



We now discuss in detail the choice of function in Eqs. (|9]) and (10). The function ijj]^ 
as defined in Eq. (10), transitions rapidly from ipi (0) = {ip^ (1) 
a width 



~h 



1 



civ 



or 



T^G 



\R\ 



0) to unity, across 
(13) 



where the value of Ai {A^) is related to the width of ipi (ipw)- Note also the existence of the 
scale 

h 

(14) 



ReV 



which is the channel midpoint where dU/dz = 0. The choice of wall function ipi rapidly 
dampens the eddy viscosity to zero near the interface, but has little effect elsewhere. Using 
this choice, together with the form 



G{s) = s{l- s)V{s) , 
we obtain the correct viscous behaviour for the velocity profile near z = 0: 

U ~ Const. + — - [1 + (s)] ds, as z ^ 0, 
Reo Jo 



(15) 



TD„2 

Const. + + O {zf 

Reo 



as 2; — )■ 0. 



(16) 



The form of G given in Eqs. ^ and (15) is chosen such that the basic velocity profile 
possesses a log layer close to, but not at the interface (wall). The G-function we use (Eq. ^) 
was derived in the paper of Biberg (2007). Our model generalizes this work by taking account 



of the dynamically important viscous sublayers. This extra detail has the added advantage 
that logarithmic singularities are no longer present in the velocity profile. Thus, for those z- 
values in the part of the domain sandwiched between the interface and the channel midpoint, 
that is, for 



1 



Re. 



< 2; < 



Rel 



ReV 
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the function G (z) has the property that 



Giz) 




+ 



(17) 



and hence, 



U ~ Const. + 



Re 



ds 



= Const. H — ^ In (z) . 

nRe 

A similar calculation near z = 1 establishes the existence of a log layer close to the upper 



wall. Finally, we establish the values of Ai^^ and n in Eq. (10). Close to the interface z = 0, 
the Reynolds stress has the form 



ttss 



Rel 



PcUi Rel A, ' 

Now a good approximation to the interface in two-phase turbulence with a large density 
contrast is in fact a solid wall (Fulgosi et al. , 2002 ). Thus, we set n = 2, the value appropriate 
for wall-bounded turbulence (Pope, 2000). In a similar manner, we fix Ai and A^ with 



reference to single-phase theory, wherein there is a one-to-one correspondence between the 
values of A^ and A^ and the additive constant B in the single-phase log law. With k = 0.4, 
the specification 



Ai,^ = e^'-^Re^.^ (18) 

corresponds to the known value B = 5.3. It is this relationship that we use throughout 
our study. We combine these modelling assumptions to obtain a velocity profile in Fig. |2} 
This velocity profile is computed for Reo = 1000, for which the corresponding superficial 
Reynolds based on the gas flow rate is approximately 12, 000. The near-interfacial viscous 
and logarithmic layers are visible in Fig. |2] (b). 

The numerical solution for the base-state profile also enables us to determine the friction 
Reynolds number and the liquid Reynolds number as a function of the control parameter 
Reo'. this is done in Fig. [3j This liquid Reynolds number is defined as 



Re, 



(19) 



Figure [3] (a) shows the dependence of the friction Reynolds number i?e=K on the control 
parameter Re^. The relationship is approximately linear. This is a consequence of the 
very small velocities in the liquid, compared with the maximal gas velocity. Thus, the gas 
layer closely resembles single-phase channel flow, and the condition U (0) <^ t^max mimics 
the zero interfacial-velocity condition in the single-phase channel. The channel midpoint 
where dU/dz = is thus approximately equal to half the gas-layer depth h^^ ^ h/2. Using 
this guess in Eq. (14) gives Re^ ~ Reo/\/2, which is close to the slope calculated in the 



figure. Figure |3] (a) shows the dependence of the liquid Reynolds number ReL on the control 
parameter Reo- The approximately quadratic relationship is a consequence of the definition 
of ReL (Eq. ([l9|), and the linear relationship between Reo and i?e^,. 
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(b) 



FIG. 3: Properties of the base state as a function of the Reynolds number Rcq. We have set 
6 = 0.1 and (m, r) = (55, 1000). (a) Dependence of the friction Reynolds number i?e* on Reo] (b) 
the Reynolds number Rbl. 



Symbol 


Numerical value, 


S.I. Units 




1.8 X 10" 


-5 


m = ni/fJ-G 


55 




PG 


1 




r = Pl/pg 


103 




dL 


10-3-10" 


-2 


6 = di/h 


0.1 




9 


9.8 




a 


0.074 





TABLE I: Table of parameter values used to estimate the wave speed. 



B. Slov^ and fast waves 

Our model also gives a way of predicting the values of (-Rcq, di) for which the critical- 
layer instability could be relevant. This mechanism depends sensitively on the shape of the 
base state, and causes a tiny wave-like perturbation at the interface to grow in time when 
{d'^Uc / dz^) ^^^^ < 0, where the critical height z^. is the root of the equation Uq (z) = c, 
and where c is the wave-propagation speed. When the critical height lies inside the viscous 
sublayer, the curvature of the mean profile is negligible, and this mechanism is not relevant. 
We obtain an estimate for the Reynolds numbers Reo that produce this regime where the 
critical layer is unimportant by solving the equation 

Rc^ Rc^ _j_ ^ ^ J- 

-Zq = — z^ — c, ^ o, 



Ren Ren 



c Re^ 

TT ^ ^1^- 20a 
Uq Reo 
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FIG. 4: An estimate of the boundary between slow and fast waves, as a function of Reynolds 
number Reo and wavenumber a. For very thin films {di = 1mm), slow waves are guaranteed at 
almost all Reynolds numbers, while for thicker films {di = 10 mm) the waves are faster for all but 
the highest Reynolds numbers. 



Equation (20a) gives a formal definition of a slow wave. Since i?e* 
layers, this definition reduces to 

5 



0(1) 



Reo/ V 2 for thin liquid 



(20b) 



We estimate the wave speed c, which we denote by Cest, by using the formula for gravity- 
capillary waves on a quiescent free surface (recall that the tilde is used to denote dimension- 
less variables): 



Cest 



gh 



r - 1 1 
- + 



a 



\/ tanh [aS). 



(21) 



We test the accuracy of this formula in a number of cases (see Sec. Ill B ): it gives an order-of- 
magnitude prediction of the wave speed. We use the values from Tab. |I] and obtain a graphical 
description of the boundary between slow and fast waves, as a function of the parameters 
{dL,S,Reo,a) (Fig. |4]). When Cest/^o ^ 1; we expect the critical-layer mechanism to be 
unimportant. This is precisely the regime of small di-valnes and high Reynolds numbers, 



which is the subject of this report. According to the classification of Boomkamp and Miesen 
(1996), the other two mechanisms of instability that exist for two-phase flow are the viscosity- 
driven instability, and the liquid internal mode. We must therefore be on the lookout for 
these instabilities in the linear stability analysis that follows. 
One final note concerning Eq. (21 ). It can be re- written as 



Cest 



Fr 1 

- + 



S 



-a 



r + la r + 1 

'ghr -11 5 
Uq r + la r + 



A/tanh {a6), 



-ctA/tanh {a6), 



(22) 
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where we have introduced the inverse Froude and inverse Weber numbers, respectively 



Fr 



g {pL - Pg) h gh 



PcUi 



r/2 



r-r 



s 



a 



PcU^h 



(23) 



By varying the inverse Froude number Fr 
accomphshed. 



a transition between slow and fast waves is 
Equation (|22|) makes the role of Fr manifest in this process. Moreover, 



it suggests the possibility of generating fast waves by fixing Fr and reducing the density 
contrast, or increasing the strength of the surface tension. We shall return to this question in 



the parameter study in Sec. IV It is important, however, to treat this analysis as preliminary, 
since we have no right to assume that Eq. (21) is valid. Indeed, a central message of § IV is 



that the wave speed must be determined, along with the growth rate, by an Orr-Sommerfeld 
type of analysis. 



C. Interfacial roughness 

So far we have been concerned with flow profiles where the interface is a perfectly smooth 
surface separating the phases. Now, we allow for surface roughness by modifying the eddy- 
viscosity law Q and Q. The work of |Lin et al.] ( |2008D gives one possible explanation for 
the generation of such roughness. This work indicates that the so-called Phillips mecha- 



nism (Phillips, 1957) may be important, whereby instantaneous turbulent pressure fluctua- 



tions give rise to a regime of linear wave growth. This is later followed by an exponential 
growth regime, which is primarily governed by the disturbances in the flow induced by the 
waves themselves. In the present context, we regard the surface roughness as a consequence 
of the gas-phase turbulence, which then acts on the interfacial waves, thereby modifying the 
growth in the wave amplitude. 

Two approaches to the modelling of the surface roughness present themselves. The first, 



and more rigorous approach, is to use an eddy- viscosity model like that of Biberg (2007). 



Such a model has the effect of shrinking the viscous sublayer near the interface. In our 
formalism, this is achieved by altering the form of the mixing length near the interface: 
before it was 

C ~ zifji (z/h) , as z — )■ 0, 

where ipi (z/h) is the damping function that operates in the near-wall region z < ^U^-Jpg') 
now, instead, we propose the behaviour 



as z — )■ 0. 



Thus, 



^ = -^[^ + i^(l-5)](l-^)^w(l-^)V(.)^, 



z/h 



d^' 



where K = £[/ {nh) is the nondimensional interfacial roughness parameter. Now the second, 
and more ad hoc approach is simply to reduce the interfacial viscous sublayer region in our 
fiat-interface model. This is accomplished by reducing the parameter Ai in the wall function 
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FIG. 5: A comparison with the work of Willmarth et al. (1987) for single- phase channel flow 
for Rcra = 2130 X 10^. Excellent agreement between the model and the experiment is obtained, 
throughout the flow domain. 



i/^i. Morland and Saffman (1993) have used this approach, and have parametrized the effects 



of roughness simply by reducing the viscous-sublayer region of the flat-interface model. 
Such a reduction then gives rise to a reduced growth rate. This is in itself a rather trivial 
observation, although it does have important implications for the linear stability analysis 
of two-phase flow: when the interfacial growth rate is reduced, there is the possibility of 
mode competition, and an internal mode can come to dominate in the stability analysis. 



We comment on this in Sec. IV C where we compare and contrast the results of a linear 
stability analysis using both roughness models. 



D. Comparison with other studies 



To validate our model for the base state, we compare it with other studies of both 
single- and two-phase flow. A further comparison with studies of single-phase flow over a 



wavy wall (Abrams and Hanratty 1985 Zilker et al. , 1976) is provided in Appendix A. We 
flrst of all characterize the single-phase version of our model. This is obtained by setting 
Uint = and by ignoring the liquid layer. We compare with the experimental work of 
Willmarth et al. (1987), for single-phase pressure-driven channel flow. In the experiment, 
the Reynolds number based on the friction velocity was 1.143 x 10^, which corresponds to 
a model Reynolds number Rcq = a/2 x 1.143 x 10^. The mean Reynolds number in the 
experiment was i?em = 2.158 x 10^ - a Reynolds number based on the mean velocity and 
channel depth. We compute Rem to be 2.13 x 10^, close to the value given in the experiment. 
A plot of the proflle is shown in Fig. [5j The model and the experimental data are in excellent 
agreement. 

To validate the two-phase version of our model, we flrst of all compare it with the work of 



Akai et al. (1980, 1981) whostudied two-phase turbulence for an air-mercury system, where 



m 



77, 



1.120 X 10' 



at room temperature. 



The liquid Reynolds number, based on the liquid-layer depth and the mean liquid velocity is 
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(a) 



(b) 



(c) 



FIG. 6: A comparison with the work of Akai et al. Subfigures (a) and (b) show the comparison 
for Rera,G = 2340; (c) shows the result for i?em,G = 3690. Here U{l/2) denotes the hquid-phase 
velocity half-way between the bottom wall and the interface. In both cases, -Rcm.L = 8040. The 
agreement is excellent in the gas phase and reasonable in the liquid phase, and is almost identical 



to the predictions given in the paper of Biberg (2007) 



set to -Rcm.L = 8040 throughout the experiments. Because the liquid is no longer laminar, we 
apply the turbulence model Eq. ^ to both layers. The gas-layer Reynolds number, based on 
the gas-layer height and mean gas velocity, varies: we study the cases where -Rcm.G = 2340 
and 3690, for which the interface is flat. We have obtained the value of 6 corresponding to 
the flow rate Rera,L = 8040 through numerical iteration^. In a similar way, we have obtained 
the value of Rcq corresponding to the flow rates Rem,G = 2340 and 3690. The results of 
this comparison are shown in Fig. [6| where excellent agreement is obtained, particularly in 
the gas phase. The agreement between the model and the experiments is as good as in the 



paper of Biberg (2007). This is not surprising, since our model is designed to replicate his 



in the log-law regions of the flow, and in the core regions. Indeed, we conclude from the 
near-exact agreement between our predictions and those of Biberg that our model inherits 
all the results he obtained from experimental comparisons. The added advantage of our 
model is that it can be continued down to the wall and interfacial zones. 

Finally, to validate the near-interface region of the model, we compare it with the DNS re- 



sults of Solbakken and Andersson (2004) for two-phase lubricated channel flow. To compare 
with their results, we take S = 1/34, m = 2, r = 1 (hence Pl = Pg = p); and 



p{h + 2dL)Ur p{h + 2dL) h + 2dL 



2p 



G 



2p 



G 



2p 



dp 



dx 



180. 



Thus, 



25 + 1' 



Re^ 
Re.j 



25 + 1 



3/2 



^ We use two parameters to specify the flow configuration: 5 and Re^. In many experiments, liquid and 
gas growth rates are used instead. The latter can be obtained from the former within the framework of 
the model through numerical iteration. 
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FIG. 7: Comparison with the results of Solbakken and Andersson (2004) for a lubricated channel. 
The broken-line curve gives the model profile across the channel, while the solid-line curve describes 
the DNS results. The latter results only extend to the channel midpoint. 



Furthermore, we take 



z 2Re, 



u- 



U_ 



h2S + r t/o V 2(5 + 1 

The results are shown in figure [7| There is excellent agreement, in particular near the 
interface. We have also compared the model predictions with the two-phase numerical 
simulations of Adjoua and Magnaudet (2009) for air and water (not shown), and have found 
similar good agreement. Thus, our model is an adequate base state for the stability analysis 
we now carry out. 



III. THE MODEL PERTURBATION EQUATIONS: DERIVATION AND 

NUMERICAL STUDIES 

A. Derivation of the model perturbation equations 

We base the dynamical equations for the interfacial motion on the Reynolds-averaged 
Navier-Stokes (RANS) equations. The turbulent velocity is decomposed into averaged and 
fluctuating parts {U, W) and («', w') respectively. The averaged velocity depends on space 
and time through the RANS equations: 

V • t/ = 0, (24b) 

where (■) denotes the averaging technique. We use these equations to model a flat-interface 
or base state of the two-phase system shown in Fig.[l] Next, we introduce a small disturbance 
that shifts the flat interface ai z = to z = r] (the dimensionless wave elevation), where 
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\v\ < 1- 

follows: 



This induces a change in the average velocity and pressure fields, denoted as 



{U, W, P) = (f/o (z) + 6u (x, z, t) , 6w (x, z, t) , Pq (x, z) + 6p (x, z, t)) 



where we denote base-state quantities by a subscript zero. Since the fiow is turbulent, and 
since the perturbations take the form of a wave, they must satisfy the RANS equations for 
a linear wave, dt = —cdx'- 



{Uo - c) -^5u + ^^Sw 
ox dz 

d 

piUo-c) —6w 

9 , d ^ 
—du + —dw 
ox oz 



The quantities 

5r = -(MV)-r(°), (5(T, 



d f d'^\ d d 

d f d 



0. 



-{w''^)-af \ 6a = 6a^-6> 



are the perturbation stresses due to the turbulence in the perturbed state, while the quan- 
tities with the zero-superscript are base-state stresses. Using the streamfunction represen- 
tation {5u,5w) = {(l)z,—4>x), and the normal-mode decomposition dx = ia, the perturbed 
RANS equations reduce to a single equation. In non-dimensional form, the gas equation is 



{Uo - c) (D^ 



a 



2\ A ^0 



dz2 



(D^ - a'^f (t)G + iama+ (D^ + a^) 6t, (25a) 



where D = d/dz, while the liquid equation is simply 
(f/o - c) (D^ - «2) 0^ 



d'U, 



dz^ 



m 



Reo 



2\2 



(D^ - a') 



(25b) 



Orr 


1907a|b 


Orszag 


1971 



Yiantsios and Higgins, 1988), with extra turbulent stresses in the gas. The problem of 
modelling the additional stresses in Eq. (25a) is considered throughout the literature. In 



this section, we use two stationary turbulent models from this literature to describe these 
stress terms. Both models give the same result. A stationary model is appropriate for slow 
waves because the dynamically important region for the instability located very close to the 
interface, on the gas side. This is precisely the region of the gas domain where the turbulence 
stationarity condition is satisfied, namely that the eddy turnover frequency U^i/ (kz) should 
greatly exceed the advection frequency a\UG (z) — c\ (Belcher and Hunt , 1993 Belcher et al. , 
T994l|Janssel^[2004| . 

The visco-elastic model: This is a stationary turbulence model wherein the perturbation 
Reynolds stresses are assumed to be proportional to the perturbation-induced turbulent 



kinetic energy (TKE). Such models have been used by Townsend (1972, 1980), and by lerley 
and Miles ( |2001 ). The model described here fits into the framework of the latter paper, 
with slight modifications: the base-state quantites are computed according to the formalism 
in Sec. [11} and the dissipation rate is taken to be linear in 6k. This last assumption is 
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not necessary, but is convenient from a mathematical point of view, since in this form. 



the dissipation rate is well-behaved at the interface, unlike the other models (lerley and 



Miles, 2001 Townsend, 1972, 1980). The perturbation- induced TKE satisfies a balance 



law wherein the advection of the kinetic energy is balanced by production, dissipation, and 
diffusive effects. The production of TKE is proportional to the stresses 6ri2 = —Sr and 
6r = —6a, the dissipation term is assumed to be linear in 6k, while in the dynamically- 
important near-interfacial region, the molecular viscosity is expected to dominate over the 



turbulent viscosity (whose effects are anyway always negligible (lerley and Miles, 2001 
Townsend, 1972 1980[ )). Thus, we have the following balance law: 



ia {Ug - c) + 



Rei 
Ren 



6k 



Re, 







a') 6k-6r,,^-ru {D' + a') 



To close the system, the visco-elastic hypothesis is invoked: 



6ri2 



-^6k 
kc 

^G 

6r - —6k = 
kG 



0, 



0. 



— zarD</)G' + ia— (26a) 
az 



(26b) 



(26c) 



The base-state stress r = — rfg is modelled as in Sec. [n| and the stress r is set equal to /cq, 
the base-state kinetic energy, consistent with the DNS results of Spalart (1988). Finally, the 
base-state turbulent kinetic energy is modelled by the equation 



kn 



PGUi 



1 Re 
C^R^l 



1^ {~z)^il-z) 



(26d) 



where C is another constant, here taken to be 0.55, which is the value appropriate for the 
logarithmic region of the mean velocity in a boundary layer. This form is chosen because it 
accurately models the viscous and logarithmic zones in single-phase flow. 

The zero- equation model: We shall compare the results of the visco-elastic study with 
the results for an eddy-viscosity model. In this formalism, the normal stresses are set to 
zero, and the shear stress is modelled as 



6t = i^T (D^ + a^) 



where 



-Giz)ij, (5) ^„(l-5) 



(27a) 



(27b) 



as in Sec. |TT} This is a rather basic model of the eddy viscosity, which does not take ac- 
count of perturbations in the eddy viscosity function itself. In particular, small changes 
in pressure will modifiy the Van Driest coefficient Ai, thus changing the viscous-sublayer 



thickness (Zilker et al. , 1976). This contribution is expected to be negligible in our small- 



amplitude analysis, and our finding that the details of the instability depend on conditions 
at the interface, and that a small modification in the extent of the viscous sublayer has 
little effect, strengthens this contention. Moreover, Belcher and co-workers have used a 
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similar model to Eqs. (27), where it was thought to capture the physics of the equilibrium 
turbulence. 



To close Eqs. (25), they are matched across the interface z = 0, where we have the 



following conditions: 

(f^L = (f>G, (28a) 

m (D^ + a^) (f)L = (D' + a^) (f)G, (28c) 

m (l)^(f)L - 3a^D(f)L) + iarRe (c - Ul) D0l + iarRe^(f)L - ^"^ff (Fr + a'^S) (t)L 
^ ' dz c — Ul^ ' 

= (D^0G - Sa^D^c) + ia^e (c - Ul) D0g + iaRe^<pG- (28d) 
The no-slip conditions are applied at 2; = —d^ and z = h: 

(Pl i-dL) = i-dL) = (pG (h) = D0G (h) = 0, 
and, where applicable, the following conditions are applied to the turbulent kinetic energy: 

Sk = z = 0, and for z = h. (29) 
Finally, at the top of the gas domain z = h, we have the no-slip conditions 

<Pg (h) = D0G {h) = 0. 



The OS equations (25) and the turbulence model reduce to an eigenvalue problem in the 
eigenvalue c. This is solved numerically according to a standard method, described and vali- 
dated elsewhere by the current authors in the context of absolute and convective instabilities 



in laminar two-phase flows (Valluri et al. , 2010). 



B. Preliminary numerical studies 

To compare the turbulence models, we carry out a stability analysis based on the values 
in Tab. [T], with di = 2.5 mm and 6 = 0.05. The inverse Froude and Weber numbers are 
computed as 

gh r — 1 , g^r — 1 



qh r — 1 , r.. r — \ 

Fr = ^—^ = (3.7809 x 10^) — ^ 

ifJ'G/pch) -f^eo Rcq 



S = L = 1-1420 X 10^ 

fil/pchRel Rel " ^^^^ 

We select a Reynolds number that produces substantial shear in the liquid, but is such that 
the liquid remains laminar. Thus, Rep = 5 x 10^, Reo = 1000, and Rcl = T\dLl iJiL ~ 460. 
According to Fig. |4| these values of (-Reo, S.di) will produce the slow- wave instability. We 
verify a posteriori that the liquid remains laminar, in the sense that the TS mode associated 
with the liquid has a negative growth rate. This growth rate can become positive, however, 
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FIG. 8: Test bed for the numerical solver. The parameters are given at the start of the sec- 



tion (IIIB). (a) Long- wave analysis of the OS model (25), without the PTS. The numerical wave 
speed at Q = agrees with the analytical formula, obtained by solving the OS model analytically 
at a = 0. At lowest order, the growth rate is zero, (b) Short-wave analysis of the OS model. The 

oo. The maximum dis- 



numerical wave speed agrees with Eq. (21) for free-surface waves, as a 
crepancy between the free-surface estimate and the true value of the wave speed is in the unstable 
part of the spectrum. The results in (b) hold over a wide parameter range: in (c) we compare 
the formula with the numerical calculations for a range of Fr- values. The agreement between the 
calculations confirms the correctness of our numerical technique. 



when mode competition takes place; this is the subject of Sec. |IVB[ Although our numerical 



method has been validated elsewhere (Naraigh and Spelt, 2010 Valluri et al. , 2010) we 



provide a further quick by comparing the zeroth-order long-wave analytical solution to the 
OS equation with our numerical method. At lowest order in a, the OS equation without 
the PTS reduces to the following matrix problem: 
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with solution 



Co — f^int + 



25'Rel {m-6-l+m5) 
Reo {6^ + AS^m + 6S'm + AmS + m') 



(31) 



Unlike in studies of laminar flow (Sahu et al. , 2007 Yiantsios and Higgins, 1988), the first 



order term is not available explicitly, since the complicated (i.e. non-polynomial) form for 
the base state precludes an explicit solution to the first-order streamfunction. Nevertheless, 



the formula (31) serves as an adequate test for our numerical scheme, as demonstrated in 
Fig. Is] (a). One further test is to check that the wave speed Cr agrees with the free-surface 



formula (21) in the limit of large a. This is shown in Fig. [8] (b). This figure also vindicates 
the use of the free-surface wave-speed formula in our estimate for the parameter range where 
the critical-layer mechanism is important. 

Having validated our numerical technique, we turn to a full-spectrum calculation. In 
Fig. [91 we obtain the growth rate using three models: the basic OS equation without the 



PTS (the so-called quasi-laminar approach), the visco-elastic model (26), and the eddy- 
Over a large range of Reynolds numbers (-Reo 



500-5000, Rei 



viscosity model (27). 

10^10^), the growth rates for the different models differ only quantitatively. 
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FIG. 9: Comparison between the models for the PTS. Sohd hne: quasi-laminar model; dashes: 
eddy-viscosity model; dots: 'visco-elastic' model. Shown is the parametric dependence of the 
growth rate on the Reynolds number for (a) Reo = 875; (b) Reo = 1000; (c) Rcq = 2000; (d) 
Reo = 5000. We have set m = 55, r = 1000, and 6 = 0.05. In each case, the difference of the 
maximum growth rate between the models is less than, or equal to 10%, which justifies the choice 
of the quasi-laminar model throughout the rest of this work. 

In particular, the differences between the quasi-laminar calculation and the eddy-viscosity 
calculation are small: the shift in the maximum growth rate upon including the PTS is less 
than 10% in the cases considered here, while the cutoff wavenumbers are virtually unchanged. 
The differences between the quasi-laminar calculation and the visco-elastic calculation are 
slightly larger. In particular, the cutoff wavenumber is shifted to a higher value in the 
Rcq = 5000 case (Fig. [9] (d)). Nevertheless, the shift in the maximum growth rate upon 
including the visco-elastic terms is not more than 10% in the cases considered in Fig. |9j The 
minor discrepancy in behaviour between the visco-elastic model and the other two models is 
due to the lack of understanding in modelling the kinetic-energy dissipation function, here 
assigned the simple linear form (Rel/Reo) 6k. Accurate modelling of this term will be the 
subject of future work. 

Our conclusion from the small differences evinced by these comparisons is that we are 
justified in considering the quasi-laminar approximation for the rest of this work. Further- 
more, we can provide a physical justification for the smallness of the contribution made by 
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the PTS in the eddy-viscosity modeL We use the analogy between the Reynolds-averaged 



Navier-Stokes equations and the equations for a laminar non-Newtonian fluid (Zou, 1998). 



In the latter case, a linear stability analysis has been performed on a two-phase stratified 
flow, where the bottom layer is a Bingham fluid (Sahu et al. , 2007). In that case, the au- 
thors found a small difference (~ 10%) between the results of the Orr-Sommerfeld analysis, 
depending on whether the perturbation non-Newtonian stresses were included. The differ- 
ence in the stability results was driven by the presence of extra terms in the perturbation 
equations for the bulk flow, and by the existence of extra terms in the interfacial conditions, 
which enhance the viscosity contrast. In our case, the additional terms in the bulk equations 
scale as nRe^:/ Rcq, which for thin layers is approximately k/\/2, and thus has a small effect. 
Moreover, in our case, the additional interfacial terms are turbulent in nature, and are thus 
damped to zero in the viscous sublayer, and vanish at z = 0. Hence, this second contribu- 
tion to the modifled growth rate is also small. Thus, in the case of equilibrium turbulence 
considered here, the effects of turbulence are felt almost entirely through the choice of base 
state. The only possibility for the effects of turbulence to enter through the PTS is when 
the critical-layer instability is present, since then rapid distortion effects are possible, and 
an interaction occurs between the turbulence and the critical layer. This will be the subject 
of future work. In the present paper however, we consider slow waves, and thus the effects 
of turbulence need be considered only in the base state. 

Next, we present the growth rate only for the quasi-laminar model in Fig. (10), where 
6 



m 



55, 



1000, 



0.05, and Rcq = 1000 {Rei ~ 460). Maximum growth occurs at 
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FIG. 10: Growth rate and wave speed for pressure-driven channel flow, with turbulent base state. 
We have set Reo = 1000 (hence Rei ~ 460), m = 55, r = 1000, and 5 = 0.05. The most dangerous 
mode is at a w 20, which equates to a wavelength i ^ l.Oti^,. 

a wavenumber a ~ 20, that is, for a wavelength i/di = 2tt/ (205) ~ 6. The wave speed 
Ct/Uq is less than unity for the dynamically important unstable waves. Thus, according to 
the analysis of Sec. [TT| neither the critical-layer mechanism, nor rapid-distortion effects, will 
be relevant. Note flnally the convexity of the small-a part of the growth rate in Fig 
The corresponding part of the growth rate is concave for shear-driven flow 



see 



10 



Miesen and 



Boersma (1995)). Having validated and compared our models, we identify the source of the 



instability and investigate its dependence on the various parameters in the system. 
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IV. LINEAR STABILITY ANALYSIS 

In this section we present detailed results of the Orr-Sommerfeld (OS) analysis, based on 



the reference values of the inverse Froude and Weber numbers already described in Eq. ( 30 ) . 
First of all, by performing an energy-decomposition, we confirm that the instability at work 
is the viscosity-contrast mechanism. This decomposition or budget is obtained from the 



RANS equations, and was introduced by Boomkamp and Miesen (1996): 



Tj (^^^5uj + Uf ■ V5uj + 5uj ■ VUf'^ = V ■ 5J^^^ - r^V • 5r^^\ (32a) 
^j^i -Sp + fidju fi {djw + dzSu) \ ^ / -(5rii + 5r22 -5ri2 \ ^^^^^ 



H {dx5w + dz5u) —5p + ^d^Sw J \ —Sri2 

V ■5uj = 0, (32c) 

which we multiply by the velocity 6uj and integrate over space. We obtain the following 
balance equation: 

KINj = REYj + Y DISS, + TURBj + INT, (32d) 

j=L,G j=L,G j=L,G j=L,G 



where 



KINj = dx J dzTjdu], (32e) 

REYj = -rj I dx I dz5uj5wj'^, (32f) 



DISSj = -'^ I dx I dz 



TURBj = 5j,G< r I dx I dz 
Additionally, 



' d f d d 

6r—6u + 6^12 iz-Su + -7—5w 



dx V dz dx 



(,32g) 
(32h) 



INT = j dx [6ul61l,zx + ^wl^l,zz]-^=q - j dx [6ug61g,zx + wg^'^g,^ 

which is decomposed into normal and tangential contributions, 

INT = NOR + TAN, 

where 

NOR = j dx [5wlSJl,zz - 5wg5'Tg,zz]^=q , 

and 

TAN = J dx [5ulSTl,zx - 5ug5J g,zx]^=q ■ 
For the quasi-laminar model under consideration throughout this paper, the term TU RB is 



set to zero. We perform the energy decomposition for the most dangerous mode in Fig. 10 
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KINg 


KINl 


REYl 


REYg 


DI55l 


DISSg 


iVOi? 


TAN 


0.18 


0.82 


2.34 


-11.90 


-4.28 


-57.42 


-2.73 


74.99 



TABLE II: Energy budget for the most dangerous mode a = 20 at Re^ = 1000 (hence Rcl ~ 460), 
m = 55, r = 1000, and 5 = 0.05. It is the TAN term that gives rise to a net positive energy, and 
thus destabihzes the interface. 




FIG. 11: (a) The phase shift between the viscous shear stress at the interface, Txz ix,z = 0) and 
the interface shape rj (x) for the most dangerous mode a = 20. The shift is smaU 6(j) = 0.22 x {it/2). 
Thus the tangential term in the energy budget is positive (destabilizing), as required by ( p3| ); (b) 
the behaviour of the viscous shear stress across the interface. 



{a = 20), and demonstrate the result in Tab. [Ilj The results in this table indicate that it is 
the TAN term that is the main source of the instability. This is the work done by tangential 
stresses on the interface, and can be written as 

TAN = f dx [{6uL - Sug) SJ^^l^, , i = 27i/a. 
Jo 

Since the kinematic condition implies the following jump condition on 6ul — Sug, 

6ul - Sug = V [U'g^i - Ul) = vU'l {m - I) , on z = 0, 
the tangential term is non-zero because m 7^ 1: 

TAN = (m - 1) — *- / dxT] (x) 6Jx, (x, z = 0) . (33) 
Reo Jo 

Thus, the viscosity contrast m > 1 induces instability, provided the interfacial shape t] {x) 
and the disturbance stress T^z {x,z = 0) possess a phase shift in the range [— |, |] (see 

Fig.m. 



So far we have chosen a set of parameters that are comparable in magnitude to those 
describing an air-water system under particular conditions. However, we wish to quantify 
the stability properties of the system in full generality (and in particular, to delineate the 
boundary between slow and fast waves), and we therefore investigate the implications of 
varying the pressure gradient, the density contrast, and the Froude and Weber numbers. 
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A. The interfacial mode 



In this section we set Rep = 5 x 10^ {Rgq = 1000). Fig. 12 shows the result of varying 
the parameter r through a range r = 10-10, 000. For large r-values, the maximum growth 
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FIG. 12: The effects of varying the density ratio r. We have set Re^ = 1000, {m,5) = (55,0.05), 
and we have taken Fr = 3.7809 x 10^ (r - 1) / Rel and S = 1.1420 x lO'^/fleg. (a) Decreasing the 
density contrast is destabilizing; (b) the dependence of Amax on r; (c) the dependence of the wave 
speed on r. Decreasing r leads to faster waves, although the unstable waves are still slow. 

rate decreases upon decreasing r. This is not surprising: a decreasing value of r implies that 
the density of the liquid approaches that of the gas, and thus the liquid has less inertia. 



The interface is then expected to be less stable. The plot in Fig. 12 (b) neatly sums up 



this dependence. However, for smaller r-values r < 10, this dependence is reversed. This is 



explained by the energy budget in Tab. Ill, where the energy decomposition at the maximum 
growth rate is shown, as a function of r. The principal source of instability is the TAN term, 
consistent with a viscosity-contrast instability. Note that the term REYl is positive too. 
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0.00 
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3.9 
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0.01 


0.00 


0.32 
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-0.11 


5.56 



TABLE III: Energy budget for the most dangerous mode as a function of r, for Reo = 10^, 6 = 0.05, 
and m = 55. In general, there are two terms contributing to the instability: one interfacial, and 
one due to effects in the liquid layer. As r is reduced, the latter term diminishes importance, thus 
removing one of the sources of instability. 



although this contribution diminishes with decreasing r, so removing a source of instability 
and thus reversing the monotone dependence of the growth rate (at large r, the growth rate 
decreases with increasing r). As r decreases further (in particular, for r = 1, 5), there is a 
destabilizing net input of energy into the perturbations from the REYg term, which implies 
that the critical-layer mechanism plays a secondary role. 

These findings raise two questions. Is the turbulence model valid at these fast wave 
speeds? Is our base-state model vahd at these low values of r? Now when the critical-layer 
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mechanism is relevant, the dynamically important region moves into the bulk of the gas flow 
(away from the interface). It is possible that this region will coincide with that part of the gas 
domain where the rapidity of advection dominates over the eddy turnover time (see Sec. III). 
Thus, in this case, rapid-distortion effects may be important. These effects may alter both 
the growth rate and the structure of the wave-induced velocity field, and will be considered 
in future work, in a parametric study similar to this one. The second question concerns the 
wall-interface equivalence in the base-state model. This assumption was used in choosing the 
exponent in the wall function. Such an equivalence is only valid for large density contrasts. 
However, this is a relatively unimportant ingredient in the model, since it determines the 
second-order term (but not the first-order term) in the Taylor expansion of the base-state 
velocity near the interface z = 0. Thus, this equivalence assumption is unlikely to affect the 
development of fast waves as r J, 1. Having verified the effect of density stratification on the 
character of the instability, consistent with our analysis in Sec. [TT], we now perform a more 
systematic analysis on the effects of the Froude and Weber numbers on the instability. 

We perform a parameter study based on the inverse Froude number Fr in Fig. 13, wherein 
Fr is varied around the reference value Fro given by Eq. (30), at fixed density ratio and 
Reynolds number. As expected, increasing Fr is stabilizing; such an increase also leads to 




(a) (b) 

FIG. 13: The effects of varying the parameter Fr at fixed density ratio. We have set Rbq = 1000, 
{m,r,d) = (55, 1000,0.05), and have taken the reference values Fro = 3.7809 x 10^ (r - 1) /Re^ 
and 5*0 = 1.1420 x W'^/Req. This corresponds to a liquid-film depth 0.0025 m and a gas-layer 
depth 0.05 m. (a) Increasing Fr is stabilizing; for a given parameter set {m,r, 5, S, Reo) there is 
a critical Froude number for which the interface is stable at all wavenumbers; (b) increasing Fr 
increases the wave speed, in particular, the wave speed is increased in the wavenumber range for 
which instability is observed. 



faster waves, and Cr/Uo > 1 in a range of wavenumbers where the interface is unstable. 
In Fig. 



14] we demonstrate the effects varying the inverse Weber number S relative to the 
reference value Sq (Eq. (30)), at fixed Reynolds number. As expected, increasing S is 
stabilizing, and leads to faster waves. However, the fast waves are stable, while the slower 
waves are unstable. This is in contrast to Fig. [14| where increasing Fr led to fast, unstable 
waves. This difference can be readily understood by recourse to the free-surface formula (22 ), 
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3.5 




(a) (b) 

FIG. 14: The effects of varying the inverse Weber number S. We have set Rcq = 1000, (m, r, 6) = 
(55, 1000,0.05), and have taken the reference values Fro = 3.7809 x 10^ (r - 1) /Rel and = 
1.1420 X / Rcq. (a) Increasing S is stabihzing; (b) increasing S increases the wave speed, 
although this effect generates fast waves only in the wavenumber region of stable interfacial waves. 



here recalled to be 

Cr / Fr i S T / , , 

— ~ \ H -a A/tanh [aa). 

The effects of varying S in this equation are felt at large wavenumbers, while the effects 
of varying Fr are more significant at smaller wavenumbers. Since the instability attains 
maximum growth at {Ijdj}) ~ 1-10, it is the variation in Fr that is prominent in the 



Fr 


S 


a 


KINg 


KINl 


REYl 


REYg 


DISSl 


DISSg 


NOR 


TAN 


Fro 


So 


20 


0.18 


0.82 


2.34 


-11.90 


-4.28 


-57.42 


-2.73 


74.99 


Fro 


WSo 


3 


1.00 


1.00 


0.00 


-15.24 


-2.72 


-65.40 


-6.78 


91.14 


Fro 


205o 


2 


1.00 


0.00 


0.00 


-33.29 


-6.00 


-155.23 


-13.22 


208.74 


O.lFro 


5*0 


10 


0.00 


1.00 


0.00 


-10.22 


-0.53 


-39.06 


-2.22 


53.04 


lOFro 


So 


30 


-3.16 


2.16 


0.00 


0.04 


-0.89 


-1.53 


-0.63 


2.00 



TABLE IV: Energy budgets for the parameter study in which S and Fr are varied, at fixed density 
ratio r. We have set Rcq = 1000, {m,r,6) = (55,1000,0.05), and we have taken the reference 
values Fro = 3.7809 x 10'^ (r - 1) /Re^ and So = 1.1420 x W'^/Rel- In all cases considered, and 
for both slow and fast waves, the instability is due to the viscosity-contrast mechanism. 



wavenumber range of instability. Thus, it is likely to be a shift in Fr, rather than 5* 
that precipitates a change in the character of the unstable interfacial waves. In Tab. 



IV 



we verify that it is the viscosity-contrast mechanism, and not the critical-layer mechanism, 
that is at work in each case of instability studied, in spite of the change in the wave speed. 
Careful parameter tuning is thus required to observe the critical-layer instability, which was 



highlighted in Fig. 12 and Tab. Ill 



Finally, we study the effect of varying Rcq in Fig. 15 The maximum growth rate and 
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the cutoff wavenumber are both shifted to higher values with increasing Rcq. Typical values 




(a) (b) 

FIG. 15: The effects of varying the Reynolds number Reo on (a) the growth rate; (b) the wave 
speed. We have set (m, r, 6) = (55, 1000, 0.05), and have taken Fr = 3.7809 x 10^ (r - 1) /i?eg and 
S = 1.1420 X 10^/i?eg. The wave speed Cr/Uo is less than unity for the unstable waves, confirming 
that these waves are in fact slow. 



of the wave speed are higher for smaller values of -Rcq, as predicted by the formula for 
free-surface waves (21). However, for unstable waves, that is, for a less than the cutoff 
wavenumber, the wave speed c^/Uq is less than unity, confirming that these are in fact slow 
waves. Two further issues arise when studying the i?eo-dependence of the stability. First, 
upon decreasing -Reo, the lower critical wavenumber shifts from ad = to some finite value 
Oci > 0. This suggests that for a given parameter set (m, r, 6, Re^S, RcqFt), there is a critical 
Reynolds number for stability. This is demonstrated in Fig. 16, where the critical Reynolds 
number is Rcqc ~ 750. In a later section, we use this result as a means of verifying our 
model against experiments, since the critical Reynolds number for the onset of instability 
is readily measured. The second issue concerns the development of a second unstable mode 
at higher Reynolds numbers, as demonstrated in Fig. 16 (b). This is the so-called internal 
mode, which we now investigate in detail. 



B. The internal mode 



We examine the properties of the second unstable mode observed in Fig. 16 (b). We first 
of all examine the energy budget at a = 25, for both unstable modes. This is shown in 
Tab. |Vj As usual, the first mode, associated with the eigenvalue branch that has interested 
us until now, derives all but a small fraction of its destabilizing energy from the TAN term, 
which we have identified as a work done by the tangential stress on the interface. This term 
is positive when m > 1, and we designate this mode the 'interfacial' mode. The second mode 
derives the majority of its destabilizing energy from this source too, although the term REYl 
is now more important. Thus, the transfer of energy from the mean flow in the liquid, to 
the perturbation flow, is important. We therefore designate this mode the 'internal' mode. 
This justification is strengthend further by examination of the streamfunction and wave 
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(a) 



(b) 



FIG. 16: The effects of varying the Reynolds number Rcq. We have set (m, r, 5) = (55, 1000, 0.05), 
and have taken Fr = 3.7809 x 10^ (r - 1) /Rel and S = 1.1420 x 10^/i?eg. Subfigure (a) shows 
the existence of a critical Reynolds number below which the interface is stable; (b) demonstrates 
the development of a second mode of instability at higher Reynolds numbers {Rcq = 2000). 
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TABLE V: Energy budget of the interfacial and internal modes at a = 25. Here we have set 
{m,r,5,Reo) = (55,1000,0.05,2000). The values of Fr and S by Eq. ([30|. Both modes enjoy 



destabilizing contributions from TAN and REYl, although this latter contribution is much larger 
for the internal mode. 



Reynolds stress function associated with these modes, shown in Fig. 17, Figure 17 (a) shows 



the streamfunction for these two modes. The streamfunction of the internal mode possesses 
a large non-zero component in the liquid, in contrast to that of the interfacial mode. This 
gives rise to significant flow in the liquid, and hence gives an important contribution to the 
transfer term REYl. The development of a larger transfer term is shown in Fig. 17, where 
we examine the wave Reynolds stress. This is the function 



T"wivc(^) = -^i I Sui{x,z)6wi{x,z) 



dU 



(0) 



dz 



-dx, 



REYl 



''wave {z) dz. 



dL 



Clearly, REYl is much larger for the internal mode, thus confirming the importance of the 
dynamics of the liquid layer for the development of this secondary instability. Moreover, the 
critical layer of the internal mode is in the liquid, a fact that has been used in the past to 
justify its designation as 'internal' (Miesen and Boersma, 1995) [U (0) = 0.92, and Cr = 0.25 
at a = 25 for the internal mode). 

The existence of a second unstable mode implies the possibility of mode competition, 
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Intetna] 
Interfacia] 



0.05 



(a) (b) 

FIG. 17: The streamfunction for the internal and interfacial modes, normahzed such that 
max |K {(j)) 1 = 1. The wavenumber is a = 25, and the other parameters are the same as those in 
the energy- budget table |V) (b) the wave Reynolds stress function for the internal and interfacial 
modes, normalized such that max r^ave = 1- The internal mode exhibits a stronger flow in the 
liquid layer, and thus gives rise to a larger wave Reynolds stress there. 




FIG. 18: Mode competition between the interfacial and internal modes. Here {r,5,ReQ) = 
(1000,0.05,4000) The parameter m takes the values 55, 20, and 5 in subfigures (a), (6), and 



(c), respectively. The wave speed corresponding to the internal mode is shown in Fig. 19 



in which the most dangerous mode changes type, from being interfacial to internal. In 
Figs. 18-[T9] we demonstrate how this competition can be achieved by decreasing m. This 
is expected to reduce the importance of the interfacial mode relative to the internal mode, 
since TAN oc m — 1. The figure does indeed confirm a change in the character of the most 
dangerous mode as m is reduced: when m is reduced from 20 to 5, the most dangerous 
mode becomes internal. Note that this crossover depends not only on m, but also on Rcq: 
we need first of all to identify a value of Reo for which the internal mode is positive, and 
then carefully select m to observe mode competition. Now a similar modal competition has 
been observed in two-phase mixing layers by Yecko et al. (2002), where again, the mode 
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FIG. 19: Companion to Fig. 18 internal-mode wave speed as a function of wavenumber for various 
m-values, at Rcq = 4000. The continuous nature of these curves confirms that the sharp changes 
in the growth rate of the internal mode are genuine. Although the wave speeds plotted here are 
positive (to facilitate easy comparison between the different m- values), the wave speed Cr — CAnt 
is negative, confirming that the critical layer for these waves is in the liquid, thus justifying the 
designation of this instability as internal. The sharp "kinks' correspond to the turing points 
or inflection points in the associated growth-rate curves. Such kinks often occur during modal 
coalescence, for example, in the Kelvin-Helmholtz instability ( Chandrasekhar , 1961), and in other 
two- phase flow scenarios ( Shapiro and Timoshin , 2005 ) . 



competition is a function of the viscosity contrast. What these examples share is the control 
of the modal competition by a parameter that requires a change in the properties of the 
two fluids under investigation. There does however, exist a situation in which the mode 
competition can be engendered by a change in the flow properties (more precisely, a change 
in the properties of the turbulence), rather than in the fluid properties. It is to this example 
that we now turn. 



C. Modelling surface roughness 



In this section we examine the effect of surface roughness on the internal and interfacial 
modes. Surface roughness is modelled by two distinct approaches. Now it is not inconsis- 
tent to examine the effects of surface roughness on wave growth: the origin of the surface 
roughness is not found in the waves we study, but rather in instantaneous pressure fluctu- 
ations at the interface that give rise to a roughened interface, where the vertical extent of 
the roughness elements is proportional to the strength of these pressure fluctuations. Such 



fluctuations appear in the work of Phillips (1957), and direct numerical simulations by Lin 



et al. (2008) indicate that these pressure fluctuations are a precursor to the exponential wave 
growth we have described here. As mentioned in Sec. |Tl| we have two distinct models for 
the surface roughness. In the first case, we use the smooth-interface model, with a reduced 
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FIG. 20: Dependence of the growth rate on the depth of the interfacial viscous sublayer, s = 
(5,2, where d=K = i^g/U^- By comparing (a) and (b), we see that decreasing the viscous- 

sublayer thickness decreases the maximum growth rate of both the internal and the interfacial 
modes. The viscous sublayer cannot be reduced further than the value implied by = 0. 
Thus, any further increase in s beyond a certain small value s ~ has no effect on the base- 
state profile, and the growth rate is thereafter unaffected by changes in s. Here we have set 
(m, r, 6, Reo) = (55, 1000, 0.05, 4000). 



viscous sublayer thickness. Such an approach has been used before (Morland and Saffman 



1993), where it was observed that the reduced viscous sublayer produces a reduced wave 



growth rate. The second model we use is a modified version of that of Biberg (2007), where 
the eddy- viscosity contains an explicit roughness parameter K = i^/ (nh), where £i is the 
mean height of the roughness elements. We compare these two approaches in this section. 

Fig. 20 shows the effect of the viscous-sublayer thickness on the stability. The growth 
rate of the interfacial and internal modes is shown for viscous-sublayer thicknesses 5d^, 2d^, 
and ci=K. These different values are obtained from the base-state model by changing the Van 
Driest coefficient Here, = vc/U^i is the wall unit based on the interfacial friction 

velocity U^^. From the figure, we see that decreasing the viscous sublayer thickness decreases 
the maximum growth rate of both modes. 

Next, we turn to the Biberg model of interfacial roughness. Fig. 21 shows the growth rate 
as a function of the roughness parameter K . As K increases, the maximum growth rate of the 
interfacial mode shrinks dramatically, while the maximum growth rate of the internal mode 
increases slightly. This change is sufficient to promote the maximum wavenumber-growth 
rate pair on the internal branch, (amax,int, Amax,int), to the status of most dangerous mode. 



This crossover occurs for K > 0.001, as shown in Fig. [18] (b). In Fig. 21, the dispersion curve 
of the internal mode possesses a local minimum near a ~ 20, similar to Fig. 18 (c). To verify 
that this is not due to a crossover between the second and third modes, we have plotted the 
three least negative modes for K = 0.005 in Fig. 22 (a). The second and third least negative 
modes are well separated and a crossover effect is thus ruled out. Figs. 22 (b) and (c) are 
plots of the wave speed for the internal and interfacial modes: the continuity of these curves 
confirms that no crossover effect is taking place. Note, however, that the wave speeds of 
the second and third most dangerous modes intersect close to the point where the internal 



mode has its local minimum. Such phenomena often occur in modal coalescence (Shapiro 
and Timoshin, 2005). Note finally that Cr — t/jnt is negative for the internal mode, which 



shows that the critical layer is in the liquid for the internal mode. 
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FIG. 21: Dependence of the growth rate on the surface-roughness parameter, for K = 0, 0.001, 0.005 
respectively. Increasing K decreases the maximum growth rate of the interfacial mode, and in- 
creases the maximum growth tate of the internal mode, to such an extent that the most dan- 
gerous mode is internal. This crossover happens for K > 0.001, as in (b). Here we have set 
(m, r, 5, Reo) = (55, 1000, 0.05, 4000). 




(a) (b) (c) 



FIG. 22: Companion to Fig. [21] (c). (a) The growth rate of the three most dangerous modes 



at Reo = 4000 and K = 0.005, being an enlarged view of the a = 15-30 region of Fig. 21 (c). 
The sharp kink or extreme point in the growth rate of the internal mode is thus genuine, and 
not the consequence of mode crossover, (b) However, by examining the wave speed, we see that 
the wave speed of the internal and the third most dangerous modes do in fact cross close to the 
corresponding extreme point of the growth-rate curves. Finally, (c) shows the wave speed of the 
interfacial mode. 



Since the results in Figs. 20 and [21] are not identical, the two models of interfacial rough- 
ness discussed here are obviously inequivalent. Which, therefore, is the correct description? 
Reducing the depth of the viscous sublayer, and thus enhancing the extent of the logarithmic 
layer, is clearly a crude model for interfacial roughness. The level of detail in the Biberg 
model is superior, and the predictions of this model for a base state with finite roughness 



agree well with experiments, as explained in his paper (Biberg, 2007). Our prejudice is thus 
towards the latter model, and we therefore expect surface roughness to modify the stability 
properties of the system through mode competition. However, this contention must ulti- 



33 





dL (mm) 


RecH 


RecH (exp) 


Cr/Uc 


Cr/Uc (exp) 


i (inches) 


£ (inches,exp) 


(1) 


1.89 


3810 


4050 


0.13 


0.08 


1.1 


0.9 


(2) 


3.54 


2650 


2760 


0.15 


0.15 


0.7 


1.2 


(3) 


4.91 


1930 


1980 




0.19 


0.9 





a 


KINg 


KINl 


REYl 


REYg 


DISSl 


DISSg 


NOR 


TAN 


5.5 


0.99 


0.01 


4.58 


1.91 


-21.32 


-31.44 


-1.55 


48.82 



TABLE VI: Comparison with the work of Cohen and Hanratty ( 1965 ). There is excellent agreement 



between the theory and the experiments. One should note however, that the experiments carry 
a margin of error of up to 20%. Thus, the agreement between the critical Reynolds numbers is 
both indicative of the correctness of our theory, and possibly a little fortuitous. The sub-table 
is a theoretical energy-budget calculation related to experiment (3). The instability is viscosity- 
induced, although there are contributions from REYi and REYg too. 



mately be confirmed by DNS, and by experiments. Although these studies are beyond the 
scope of the present work, we are able to test the predictions of the flat-interface model 
against experiments, which we do in the next section. 



V. COMPARISON WITH OTHER WORK 



In this section, we compare our results with some of the experimental data from the 
literature, in particular the work of Cohen and Hanratty (1965), and Craik (1966). We also 



compare our findings with a model that is frequently used in practical applications to predict 
flow-regime transitions, namely the viscous Kelvin-Helmholtz theory. To do this, we refer to 
Fig. 16 (a), which highlights the importance of the Reynolds number in the stability analysis. 



In that figure, the Reynolds number is varied and the other parameters are held fixed. For 
sufficiently large values of i?eo, the dispersion curve associated with the interfacial mode is 
paraboloidal, with critical wavenumbers at ac,i = 0, and q;c,u > 0. As the Reynolds number 
decreases, the the growth rate develops an intermediate critical wavenumber a;c,uo, where 
< ttcuo < c^c.u- At a critical i?eo-value, the maximum growth rate is zero, and thus the 
intermediate critical wavenumber is simultaneously a maximum, and a zero, of the function 
Ar {a). Finally, below this critical i?eo- value, the growth rate is negative everywhere. The 
experimental works we reference involve a similar path through paramter space. 



A. Comparison with experiments 



Cohen and Hanratty ( 1965 ) report critical Reynolds numbers for millimetre-thick liquid 



films. They observe the development of two-dimensional waves above a critical Reynolds 
number. They call these waves 'fast', in the sense they move at a velocity that exceeds 
the interfacial velocity. These waves are, however, in our classification, 'slow' (or on the 
boundary beween 'slow' and 'fast'), since the theoretical values computed are Cr/f/o < 1, 
and thus the viscosity-contrast instability is expected. We show a comparison between 
the theoretical predictions of our model and the measurements of Cohen and Hanratty 
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in Tab. VI Our estimates for the critical Reynolds number Rcch = PchUc/fJ'G are in 
close agreement with the experimental values. We are mindful, however, that the margin 
of error stated in the experiment is between 10% and 20%. We have also compared our 




(a) 



1.2 
1 

^ 0.8 
0.6 
0.4 
0.2 


-0.2 







1 \ 





-0.2 



-0.1 



0.1 



0.2 



(b) 



FIG. 23: Theoretical calculation based on the parameters in experiment (5) of Cohen and Hanratty. 
(a) The streamfunction; (b) the wave Reynolds stress. In (b) we see a contribution to the instability 
from transfer terms in the liquid and the gas, although the viscosity-contrast across the interface 
gives the most important contribution to the energy of instability. 



theoretical model with the measurements of the critical wavelength and wave speed. There 
is excellent agreement between the theoretical and experimental values for the wave speed. 
The spread in values of the critical wavelength is larger, although this is acceptable, in 
view of the large error attached to the experimental measurements. The energy budget in 



Tab. VI is based on a theoretical calculation, with parameters taken from experiment (3). 
The corresponding streamfunction and the wave Reynolds stress function are presented in 
Fi g. [23j T he instability is confirmed to be due to the viscosity-contrast mechanism. 



Craik ( 1966 ) performs a similar experiment with liquid films thinner than those found in 
Cohen and Hanratty. He reports critical conditions for the generation of unstable waves. The 
trend in the data agrees with that in the theoretical calculations, although the quantitative 
agreement is poor. Craik explains that waves are observed for film thickness below that 
quoted in experiment (1.1), although the uniform thin film of liquid is difficult to maintain 
under these conditions. It is possible that the thinness of the film inhibits precision in the 
measurement at film thickness above this lower bound too. Craik also explains that accurate 
measurements of wave speed were difficult owing to the long wavelengths of the observed 
waves (compared to the channel length). These are sources of error that explain why there 
is only qualitative agreement between the theoretical and experimental data. 

We also perform theoretical calculations based on the parameters in experiments (1.5) 
and (2.3), to examine that character of the unstable waves. We provide the energy budgets 



associated with these calculations in Tab. |VIlT| The instability is interfacial: the contribution 
from REYl and REYq present in the Cohen data are absent here. This makes sense: 
REYl should be negligible because the liquid layer is so thin and thus (pL cannot contribute 
meaningfully to the dynamics, while REYq is unimportant because the waves are slow (cr/f/o 
is O(10~^) or O(10~^) for the film thicknesses and Reynolds numbers considered). 
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TABLE VII: Comparison with Table 1 (p. 375) and Table 2 (p. 378) in the work of |Craik| pi)66 ). 
There is good agreement between the theory and the experiment in the first case, and only very 
rough agreement in the second case. As explained in the experimental paper, a sharp transition 
to wavy flows was not observed in this second case, which explains these quantitative differences. 
We do not refer to experiments (2.4)-(2.5), wherein our model predicts laminar gas flow. 
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TABLE VIII: Theoretical energy-budget calculations related to experiments (1.5) and (2.1) of 
Craik. The instability is viscosity-induced, and there are no other contributions to the instability, 
unlike in the Cohen data. 



B. Comparison with Viscous Kelvin Helmholtz theory 



In this section, we compare our predictions with those obtained using viscous Kelvin- 
Helmholtz theory (Barnea, 1991). This is a simplified theory for the interfacial instability 
of two-phase turbulent flow, and takes account of turbulence in either or both phases. It is 
commonly used in one-dimensional models for large-scale stratified and slug-flow predictions. 
The velocity field enters only through the liquid- and gas-average values, ul and ug, while 
the cross-sectional area fractions = (11/ {di + (ic) and ec = ^g/ (c^l + c^g) also play a 
role. Then, the complex frequency uj is obtainable from a quadratic equation: 



o;^ — 2 {xqu — xii) CO + {x2a^ — x^a^ — ixio) = 0, 
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FIG. 24: Comparison with the viscous Kelvin-Helmholtz model: plots of the liquid film thickness 
against the Reynolds number RecmsLx = PGUcmaxdc/ IJ-G- The viscous Kelvin-Helmholtz model 
overpredicts the stability boundary by an order of magnitude. Better agreement is obtained be- 
tween the data of Craik for the 6 in channel and our theoretical model, although the margin of 
error in the experimental data is large. 
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tg 

n 



h- 

fG 



I/gPgUg, 

\fipG {ug-Ul) \ug 
(^L + (^G 



Cr. . 



Dg 



4dL, 



2d 



G, 



The coefficients Cg and Cl both take the value 0.046 for turbulent flow and 16 for laminar 
flow, while Ul and ug both take the value 0.2 for turbulent flow, and 1.0 for laminar flow. 
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Finally, the interfacial friction factor /i is assumed to be constant and equal to 0.0142 (see 



Barnea (1991)). We plot the stability boundary predicted by this theory in Fig. |24| and 
compare the results with the Craik data for the 6 in. channel (Tab. |VII ), and with a curve 
fit based on a number of points obtained from our theoretical calculations. The viscous 
Kelvin-Helmholtz model overpredicts the critical Reynolds number compared with both the 
data of Craik and our theoretical model by an order of magnitude (Fig. 24), which casts 
doubt on the usefulness of such a depth-averaged model. Our model gives better agreement 
with the data of Craik, although we are mindful of the uncertainty in these experimental 
data. Nevertheless, both our theoretical calculations and the experimental data demonstrate 
the unstable-stable-unstable transition that arises when the film depth is increased, holding 
the Reynolds number fixed. This is the statement that our theoretical curve in Fig. [24] is 



non-monotonic (the non-monotonicity in the model curve is masked somewhat by the large 
scale necessary to show the VKH results in the same figure). In conclusion, the qualitative 
agreement obtained here, together with the good agreement obtained relative to the Cohen 
data, inspires confidence in our model, while the poor agreement between our data and the 
predictions of the viscous Kelvin-Helmholtz theory calls into question the validity of this 
depth-averaged model, at least for the kind of thin-film waves studied here. 



VI. CONCLUSIONS 



In this paper, we have investigated the stability of an interface separating a thin laminar 
liquid layer from a turbulent gas in a channel. We have approached this problem in two 



steps: first by generalizing the model of Biberg (2007) to describe the interfacial and wall 



zones in pressure-driven two-phase channel flow, and second by studying the linear stability 
of this state using an Orr-Sommerfeld analysis of the Reynolds-averaged Navier-Stokes 
equations. Our model has enabled us to investigate the stability of the interface as a function 
of various parameters. In general, the interface becomes unstable due to a mismatch between 
the viscosities in the liquid and the gas. In this work, we have taken into account the 
perturbation turbulent stresses (PTS) using two distinct models: in both cases, these stresses 
have only a quantitative effect on the stability results for the thin liquid layers; for the eddy- 
viscosity model of the PTS, this effect is particularly small. These stress contributions are 
therefore ignored throughout the work. We have provided an explanation for this null result 
using an analogy with non-Newtonian fluids. This work builds upon previous work in the 
field (notably that of Miesen and co-workers (Boomkamp and Miesen 1996 Boomkamp 



et al.\ 1997 Miesen and Boersma, 1995), and Kuru et al. ( 1995[ )) by developing an accurate 



base-state model, validated against numerous experiments and direct-numerical simulations, 
and by accounting for the PTS. These efforts result in excellent agreement with the relevant 



experiments (Cohen and Hanratty 1965; Craik 1966). 



Our linear stability analyses evince a definition of slow and fast waves. A slow wave is 
one for which the instability is driven by the viscosity contrast across the interface; a fast 
wave derives its energy of instability from the critical layer. The phase speed q of a slow 
waves satisfies Cj./Uo < 1, where PgUq = h\dp/dx\. We have carried out a parameter study 
to find ways of controlling the wave speed. The inverse Froude, inverse Weber, and density 
numbers control the wave speed, as suggested by the formula for free-surface waves in a 
liquid layer. However, in all cases considered, the waves are slow. For certain values of 
the triple {r,Fr,S), there is a critical-layer contribution to the instability, although this is 
marginal. 
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For certain parameter values, we also observe a positive growth rate for the so-called 
internal mode, associated with instability that is due both to the shear content of the 
liquid and to the tangential stress at the interface. By a judicious choice of parameters (in 
particular, for small viscosity ratios), the internal mode can be made to dominate over the 
interfacial mode, and gives rise to mode competition. To engineer mode competition in this 
way, it is necessary to modify the properties of the two fluids. However, by increasing the 
level of turbulence, the fiat interface roughens, and this also has the effect of diminishing the 
interfacial mode relative to the internal mode. Morland and Saffman (1993) have explained 
previously how roughness can reduce the growth rate of the wave; here we go further and 
conjecture that this mechanism can, in addition, engender mode competition. However, this 
result is rather conjectural, and we simply mention it as a route for future experiments. 



Acknowledgements 

This work has been undertaken within the Joint Project on Transient Multiphase Flows 
and Flow Assurance. The Authors wish to acknowledge the contributions made to this 
project by the UK Engineering and Physical Sciences Research Council (EPSRC) and 
the following: - Advantica; BP Exploration; CD-adapco; Chevron; ConocoPhillips; ENI; 
ExxonMobil; FEESA; IFP; Institutt for Energiteknikk; PDVSA (INTEVEP); Petrobras; 
PETRONAS; Scandpower PT; Shell; SINTEF; StatoilHydro and TOTAL. The Authors 
wish to express their sincere gratitude for this support. 

APPENDIX A 



For further validation of the base state discussed in Sec. |TT| we compare our turbulence 
modelling with experimental data for flow past a wavy wall, obtained from the papers of 



Zilker et al. (1976), and Abrams and Hanratty (1985). The curvilinear coordinates necessary 



for this work were previously introduced by Benjamin (1958) f : 

1] = z — a$. 
If the streamfunction has the form 



(A-1) 



Uq (s) ds + aF (rj) e 



(where Uq is the single-phase version of the base state in Eq. (12)), then the momentum- 
balance equation for F is 



ta[{d'^-a')F{v)-U;;{v)F{v)]+C 
where C is the curvature-related term 



{dl-a'YF{v)+n, 



Rer ^ " 



C = 2ta'Ul,ir])Uoiv) e""" + 



1 



-Ren 



-arj 



(A-2a) 



(A-2b) 



t We thank S. Kalliadasis and D. Tseluiko for suggesting the appHcation of this coordinate system to the 
problem. 
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and TZ is the Reynolds-stress term: 



7^ = + a^) [F" (t]) + a^F (r/) + 2aU^ (r/) e""'' - 2a^Uo (v) e~"^] 

+ 2e-"^ [«7^' (r^) - aV^ (r^)] , ro (v) = (r?) f/^ (v) ■ (A-2c) 

The function fix (v) is the eddy viscosity. It is set to zero in the quasi-laminar case, and 
assigned the form of Eq. ^ if the perturbation turbulent stresses are considered. We solve 
Eqs. (A-2) subject to the boundary conditions F = F' = on ri = and on rj = 1, which 



are no-slip conditions on the perturbation F. Although the no-slip condition on the upper 
boundary is at z = 1, not rj = 1, these planes are close to one another: the physical boundary 
z = 1 corresponds to an 77- value 1 — ae~"e*"^, which is close to unity for large a- values. Thus, 
for simplicity, we impose a boundary condition at // = 1. 

The solution of the boundary- value problem facilitates a comparison with experimental 
data. In this comparison, we use the quasi-laminar assumption: the eddy-viscosity terms TZ 
are set to zero, and turbulence enters only through the shape of the base state Uo. To make 



an accurate comparison between the experiments and Eqs. (A-2) we study the shear stress 
at the interface: 



To = RCq (0^^ - (l)xx)r,=0 = 



[F" (0) + a^F (0) + 2a (Rel/Reo)] . (A-3) 



We also study the phase shift between this stress function and the wave surface a3ft [e*"^] . 
We examine the situation described by Fig. (5) in the work of Zilker et al. (1976), for which 



a/H = 0.003, aH = 13.3, Re^ = 2270, 

where H = 5.08 cm is the channel depth. We also look at Fig. (4) in the work of Abrams 
and Hanratty, where 

a/H = 0.007, aH = 27r, Re^ = 1110, 

where H is the same is in the Zilker experiment. A comparison between theory and experi- 
ment is shown in Fig. A-1 where good agreement is obtained. 

Note that in several crucial respects, the stability calculations performed in the main 
part of paper and the wavy- wall calculations performed in this Appendix are different. 
In the case of small-amplitude waves on an interface, the growth of waves is understood 
through a linear stability analysis. Thus, the amplitude of the initial interfacial disturbance 
is assumed to be small, and functions as a small parameter in a linear stability analysis. In 
contrast, the amplitude of the wavy wall is not infinitesimally small, and it is the finiteness 



of this amplitude that gives rise to the curvature terms in the equations (A-2), which in 



turn affects the distribution of stress at the interface. Furthermore, the wave speed in the 
wavy-wall calculation is a known parameter. In the linear stability analysis, it is determined 
as the solution of an eigenvalue problem. It might be tempting to guess the wave speed 



in the case of interfacial waves by recourse to the free-surface formula (21) but our linear 



stability analysis shows that the a-range of maximal wave growth is precisely that range 
where this formula is least reliable (See Fig. [s]). Thus, a key difference between these two 
calculations is that in the wavy-call case, c is a parameter; in the wavy-interface case, it must 
be determined from other parameters. One similarity between the calculations is the shape 
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FIG. A-1: (a) Comparison of the theoretical base state with the work of Zilker et aL (1976). 
We plot the total shear stress at the interface, [Rcq^Uq (0) + tq (x)] / [Rcq^Uq (0)] , and com- 
pare the theoretical curve with the data from the experiment. Reasonable agrement is ob- 
tained for the amplitude. Excellent agreement is obtained for the phase shift of the shear 
stress relative to the wavy wall: the data predicts a phase shift of approximately 50°, while 
our model predicts a phase shift 52.6°. The wavy wall is shown in the figure for compari- 
son, albeit with an exaggerated amplitude. The parameter a /Re* has the value 0.0059, while 
aa = 0.04; (b) Theoretical curves for the shape of the streamfunction and the stress distribution 
r = Re'^ [F" (??) + a^F (rj) + 2af7^ (r?) e""'' - 2a'^Uo (r?) e""'?] ; (c) Comparison of the theoretical 
base state with the work of Abrams and Hanratty ( |1985 ). The parameter a/ Re* has the value 
27r/1110 = 0.0057, while aa = 0.04; (d) Theoretical curves for the shape of the streamfunction and 
the stress distribution r. The streamfunction vanishes slowly as z ^ H , in comparison with (b). 



of the streamfunction. In the two-phase calculation, the shape of the streamfunction in the 



upper layer is similar to that obtained from the wavy-wall calculation (Figs. A-1 (b) and A- 
[l] (d)). This is a consequence of the boundary conditions, which impose severe limitations 
on the shape of the streamfunction. In conclusion, the wavy-wall calculation, because it 
assumes that the wave speed is a parameter, is an incomplete model for two-phase wavy 
interface calculations. It is, however, a testbed for verifying the turbulence modelling of the 
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